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Exact caicuiations are performed on the two-dimensional strongiy interacting, unpoiarized, uni¬ 
form Fermi gas with a zero-range attractive interaction. Two auxiiiary-fieid approaches are empioyed 
which acceierate the sampling of imaginary-time paths using BCS trial wave functions and a force 
bias technique. Their combination enables calculations on large enough lattices to reliably compute 
ground-state properties in the thermodynamic limit. A new equation of state is obtained, with a 
parametrization provided, which can serve as a benchmark and allow accurate comparisons with 
experiments. The pressure, contact parameter, and condensate fraction are determined systemati¬ 
cally vs. kpa. The momentum distribution, pairing correlation, and the structure of the pair wave 
function are computed. The use of force bias to accelerate the Metropolis sampling of auxiliary-fields 
in determinantal approaches is discussed. 
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Exact results on fundamental models are uncommon, 
especially for strongly interacting fermion systems. In 
the rare cases where they exist (for example in one¬ 
dimensional models by Bethe ansatz or density matrix 
renormalization group [11,i) , they have invariably played 
an integral role in bringing about physical insights, ad¬ 
vancing our understanding, and serving as benchmarks 
for the development of new theoretical and computa¬ 
tional approaches. 

The Fermi gas with a zero-range attractive interac¬ 
tions is a model for strongly interacting fermions which 
has generated a great deal of research activities Si- 
The model is of interest in both condensed matter and 
nuclear physics. As a model it is rather unique in that, 
thanks to advances in experimental techniques using ul¬ 
tracold atoms, it can be realized in a laboratory with 
great precision and control i i . 

In three-dimensions (3D) the interplay between experi¬ 
ment, theory and computation has lead to ramd advances 
SI. An example is seen in the evolution of the de¬ 
termination of the so-called Bertsch parameter at unitar- 
ity. Quantitative comparisons have allowed validation of 
our understanding and provided an impetus for develop¬ 
ments of both experimental and theoretical techniques. 
The remarkable level of agreement achieved recently be¬ 
tween calculation Q and experiment Q demonstrates the 
tremendous progress towards precise understanding and 
control of strongly correlated quantum matter. 

The two-dimensional (2^Fermi gas has attracted con¬ 
siderable recent interest [lll - ll^ . especially with its ex¬ 
perimental realization using highly anisotropic trapping 
potentials [l^. In 2D a bound state always exists, and 
the BCS-BEC cross-over offers rich possibilities between 
the interplay of inter-particle spacing (density) and in¬ 
teraction strength, where effects beyond the mean-field 
description will be more pronounced than in 3D. Interest 
in this model is further enhanced by the 2D nature of 
many of the most interesting and complex materials, in¬ 
cluding high-Tc cuprate superconductors and topological 
superconductors |2l| . 


In this paper, we obtain exact numerical results on the 
ground state of the strongly interacting 2D spin-balanced 
uniform Fermi gas. To date the most accurate numerical 
results on the 2D system have mainly come from diffusion 
Monte Carlo (DMC) simulations [1^. These calculations, 
however, involve the fixed-node approximation 
and lead to systematic errors which are difficult to esti¬ 
mate; furthermore, some of the correlation functions that 
are central to the physics of these systems are not read¬ 
ily available from DMC. Here, we employ two auxiliary- 
field quantum Monte Carlo (AFQMC) approaches: one 
based on the branching random walk method used in the 
3D study in Ref. Q, and the other a novel approach in 
the Metropolis path-integral framework which dramat¬ 
ically improves efficiency. Their combination allows us 
to calculate the thermodynamics and pairing properties 
exactly in the entire range of interaction strengths. 

Our calculations are performed on periodic lattices. 
We use supercells of up to 3,000 sites, containing about 
120 particles, with projection length in imaginary time of 
/3 > 50 (in unites of 1/Ep). For each lattice and Hamil¬ 
tonian parameters, the calculation is numerically exact, 
with only statistical uncertainties which are fully con¬ 
trolled. Systematic extrapolations are then carried out 
to reach the thermodynamic limit (TL). 

As the interaction in cold atoms is short-ranged com¬ 
pared to the inter-particle spacing, the uniform 2D Fermi 
gas can be modeled by a lattice Hamiltonian 

Ms 

H = t'^ eicCk^Cko- + U'^ (1) 

k,cr i 

with A/'s = sites and t = /i^/(2mA^), where A is the 
lattice parameter. Only the low energy behavior of Sk 
will be relevant, and we have used both the Hubbard 
dispersion = 4 — 2(cos kx + cos ky) and the quadratic 
dispersion = kx + ky. In this form, the momentum 
kx (or ky) is defined on the lattice, with units 27r/L, and 
kx G [—7r,7r). The on-site interaction is attractive and is 
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given by [IJ 

^ =__ ( 2 ) 

t In(kra) — In(C-yn) ’ 

which is tuned, for each lattice density n = N/Ms and 
Fermi momentum kp = ^phmjIS.^ to produce the desired 
2D scattering length a, defined as the position of the node 
of the zero-energy s-wave solution of the two-body prob¬ 
lem. The constant C in Eq. ([2]) depends on the dispersion 
relation: = 0.49758 and C« = 0.80261. 

We employ two AFQMC methods to study this model: 
the branching random walk approach, and an acceler¬ 
ated Metropolis approach with a force bias. In the 
first we project the ground-state wave function by 
importance-sampled random walks in Slater determinant 
space [ 2 ^[ 2 ^. A BCS wave function, taken from the solu¬ 
tion of the gap equation for the same discretized Hamilto¬ 
nian, is chosen as the trial wave function, and the mixed 
estimator Si! is used to calculate the ground-state en¬ 
ergy. The BCS trial wave function shortens the conver¬ 
gence time in the imaginary-time projection, and greatly 
reduces the Monte Carlo statistical fluctuations, as illus¬ 
trated in the 3D case Q. 

Our second approach is based on the ground-state 
path integral form of AFQMC, but introduces several 
advances, including accelerated sampling (described in 
more detail in Appendixby a dynamic force bias [^ . 
which enables global moves of fields on a time slice with 
acceptance ratio of over 90%, and control of the Monte 
Carlo variance [d^. Its main advantage over the the 
open-ended branching random walk approach is the ease 
with which any observables can be computed, and we use 
it to compute the momentum distribution and correlation 
functions. (Since there is no sign problem here, no con¬ 
straint is needed, which is the primary motivation for us¬ 
ing the open-ended branching random walk form.) With 
this approach, our calculations typically have /3 ^ 320 
or larger (in units of discretized with over 12,800 
time-slices. 

These technical advances result in orders of magnitude 
improvement in sampling efficiency, which makes it pos¬ 
sible to achieve the high numerical accuracy presented in 
this work. In both approaches, the computational cost 
scales as ~ The linear scaling with Afg is im¬ 

portant, as it enables calculations on large lattice sizes. 
To approach the TL, we first extrapolate calculations to 
the continuum limit by taking Afg —^ 00 while holding N 
fixed. The number of particles, N, is then increased until 
convergence is reached within our statistical accuracy, as 
illustrated next. 

Figure [T] displays the calculated equation of state 
(EOS), in units of the Eermi gas energy Epc = as a 
function of the interacting strength, x = In(fcp’a). A table 
of the AFQMC data can be found in Appendix [C] The 
top panel illustrates the convergence to the TL, where 
AFQMC energies are shown for fixed N. At each x, the 
energy has been extrapolated to the continuum limit, us¬ 
ing a 4th-order polynomial in 1/L. In the more strongly 



FIG. 1: (Color online) Calculated equation of state. The top 
panel shows the energy, relative to the final AFQMC results, 
for finite number of particles, N. Also shown are the DMC 
results of Ref. [T^, which are variational. Note the small scale 
of the vertical axis. The bottom panel shows the AFQMC 
(and DMC) results at the TL, relative to the BCS result. A 
fit has been performed on the AFQMC results for the EOS. 
The result is given in Eqs. (lip and shown as the solid line. 
The inset in panel (b) compares the calculated pressure from 
AFQMC (solid line) and DMC (dashed, taken from Ref. [^ 1 
with experiment (points) in the crossover region. 


interacting cases, we take advantage of the fact that 
and produce energies which converge to a common 
limit from opposite directions and perform both sets of 
calculations to reduce the uncertainty in the extrapola¬ 
tion. In the opposite regime, energies from the quadratic 
dispersion shows less dependence on L and they are used 
alone. We illustrate the extrapolation procedure in Ap- 
pendix|B] The error bar of each symbol, barely noticeable 
in the graph, combines the QMC statistical error (negli¬ 
gible) at each L and a conservative estimate of the un¬ 
certainty from the extrapolation, which typically involves 
half a dozen or more data points from each dispersion re¬ 
lation, with L ranging from ^ 15 to 45 (and larger if 
necessary). 

The results for different values of N show that con¬ 
vergence is reached to within our statistical accuracy by 
N ~ 100 [ 2 ^. This is consistent with DMC results [l^ 
which observed no significant change between N of 26 
and 98. The DMC results provide the current best esti¬ 
mate of the EOS and are included in Eig. I. We see that 
the error from the fixed-node approximation is largest 
in the crossover region, at intermediate values of x. The 
maximum error is about 10% of the “correlation energy”, 
the difference between the BCS and exact energies. 

In addition to serving as a benchmark for theory, the 
new EOS can provide validation for experiments. Exper- 
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iments are fast developing; in 3D remarkable precision 
0 was reached in the measurement of the Bertsch pa¬ 
rameter (with uncertainties only slightly larger than our 
symbol size in the top panel of Fig. 1). In the inset in 
the bottom panel, we show a comparison of the calcu¬ 
lated pressure with the latest experiment in 2D [2^. In 
the crossover regime, better agreement with experiment 
is seen with the new result than with DMC. There may 
be other factors contributing to the discrepancy between 
experiment and theory [^, . We leave more detailed 

comparisons of our results and experiment to a future 
publication. 

We parametrize the computed EOS by Ec = Eqmc — 
Ebcs [note that Ebcs/Efg is related to the two-body 
binding energy by I — cb/{‘^Efg), and is given by I — 
85 - 2 ( 7 - 1 - 2 :) -^^^here 7 = 0.57721 is Euler’s constant]: 

^ f fix), x< 0.2664; 

= < fix), 0.2664 < a; < 4.3058 ; 

[ .fix), X > 4.3058. 

The intermediate region is fitted with a 7th-order poly¬ 
nomial 

7 

fix) = '^aix\ (3) 

In the BCS region, the form is based on perturbative 
results ( 33 . [ 3 ^ 
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FIG. 2: (Color online) The contact parameter C. The main 
figure shows the result of C (relative to the BCS result) ob¬ 
tained from Eq. ©• The statistical uncertainty is smaller 
than the line thickness. DMC [3 and BCS results are also 
shown for comparison. The inset shows n(k)fc'^ vs fc = |k| at 
X = 0.5. The horizontal lines give the C values from DMC, 
AFQMC and BCS (top to bottom), indicated by the arrows 
in the main figure. The u(k) data are from two systems, with 
L = 45 (circles) and 51 (squares), respectively, and N = 58. 
Results are plotted for k along both the horizontal (solid sym¬ 
bols) and diagonal (open) directions. 


nx) 


4 ^ 

if. 

i=2 


( 4 ) 


while in the BEC regime a dimer form is used 


fix) = -1 + 


0.5 

'Y 


HX) Cl eLo<{^^x)^ 

X X 


( 5 ) 


where X = cq — 2x with cq = 3.703 from the dimer 
scattering length ~ 0.557a given by few-body calcula¬ 
tions 0 : and Cl = ln( 7 r) -|- 27 -|- 0.5. The parameters in 
Eqs. 0 and 0 are determined by continuity conditions 
(value and first two derivatives) from Eq. 0. The pa¬ 
rameters and the locations of the transition between dif¬ 
ferent regions are then varied in a small range to further 
minimize the variance of the overall fit with the QMC 
data. The final parameters are listed in Table mill. 

The contact is important to the physics of di¬ 

lute gase s, and can potentially be measured experimen¬ 
tally (l^. With the functional form of the EOS, it is 
straightforward to determine the contact: 

C_ ^ 1 djE/Epa) 

kp 4 dx 

The result is shown in Eig. [21 An alternative approach 
to obtain the contact parameter is from the tail of the 


momentum distribution [H, n{k)f —>■ C at large 
k. This provides an internal check on the consistency 
and accuracy of the calculation. As illustrated in the in¬ 
set, a clear plateau is present before edge effects start 
to manifest as k approaches the cut-off value, giving 
a C value in excellent agreement with that from the 
EOS. (The full momentum distribution n(k) is shown 
in Fig. [3] for three representive interaction strengths.) 
The pressure and the chemical potential can be obtained 
from simple combinations of the energy and contact: 
P/Pfg = 2Clkjp + E/Eyg^ which was applied in the 
inset in Fig. [H and ^/^fg = C/kp + E/E-pG- 

We next quantify how the pairing properties evolve as 
a function of interaction strength. The zero-momentum 
pairing matrix (of dimension Ms x A4), 

Alkk' = — ^kk'(Cjj..|.Ck-|')(c[_jj.j^C_k4,) , (7) 

is computed in the many-body ground state, where the 
pair creation operator A\^ = cjj..^c[Lk 4 ,- We associate 0 
the leading eigenstate with the pair wave function in k- 
space, ((>-|-|(k). This is shown in Fig. [3] for three charac¬ 
teristic interaction strengths. The inset shows the cor¬ 
responding real-space structures, ip-fiir), obtained from 
the Fourier transform of (/|-|■ 4 .(k). In the BEC regime, 
the momentum distribution is very broad, the pair wave 
function involves many k-values, and the pairs are tightly 
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TABLE I: Final parameter values (full digits in Appendix[D]) in the parametrization [Eqs. ([3][5])] of the exact EOS from QMC. 


i 

0 

1 

2 

3 

4 

5 

6 

7 

a\ 

ai 

al 

-11.8041 

-0.81984 

14.6755 

0.12733 

-4.85508 

0.06851 

-0.06085 

-0.01451 

0.36401 

-0.00919 

-0.61531 

0.00419 

-0.00064 

3.4312x10"® 



FIG. 3: (Color online) Momentum distribution and pair wave functions in three regimes of interaction strengths, x = In(afcir). 
In each panel, the vertical tick labels on the left are for n(k) and those on the right are for (().f^4_(k), both plotted vs. k (in units 
of kp). Note the different scales between the three panels. The inset shows the real-space wave function r in a 3D 

plot. The lattice has Ms = 2025 sites, with density n = 0.0286. 


bound like a molecule, as seen in (a). In the BCS regime 
in (c), on the other hand, modifications to the non¬ 
interacting n(k) are limited to near the Fermi surface, 
with a small number of k-vectors in its vicinity partic¬ 
ipating in pairing. The pair wave function is sharply 
peaked near the Fermi surface, and becomes very ex¬ 
tended in real space. (Residual finite-size effect can be 
seen in this case in the second ring of which is af¬ 

fected by the shape of the supercell.) As kpa is increased, 
the systems crosses over from (a) to (c) via the strongly 
interacting regime represented in (b). Beyond the central 
peak, the wave function ^/^'|- 4 ^(r) in (b) contains significant 
radial oscillations, with multiple circular nodes. 

The condensate fraction is given by the largest eigen¬ 
value of Mkk' divided by N/2. The results are shown 
in Fig. 0] as a function of interaction. At the mean- 
field BCS level Mkk' = and there is only 

one non-zero eigenvalue (equal to X]kl(^k')P)- In the 
many-body ground state, additional depletion of the con¬ 
densate is present from scattering into zero-momentum 
pairs distinct from (()-|. 4 ,(k). The BCS condensate frac¬ 
tion and pair wave functions are in reasonable agreement 
with exact results down to In(afcF) 3. For stronger 
interactions, the BCS condensate fraction grows signifi¬ 
cantly faster. At In(kpa) —I, it predicts an essentially 
100% condensate as opposed to only 80% from the exact 
result. In this regime, Bogoliubov theory of a Bose gas 
[4l| with the dimer scattering length above gives results 
consistent with the QMC data. The largest deviation 
between BCS and exact results occurs in the crossover 
region, near ln(aA:F) ~ 0.5, where the momentum distri¬ 
butions and pair wave functions also exhibit the largest 
differences. 

We also calculate the real-space on-site pairing corre¬ 


lation function: 


C'(r) = (Cot4tCr4Crt) , (8) 


where the reference point 0 and all r values related by 
translational symmetry can be averaged over. The results 
are shown as a function of r = |r| in the inset in Fig. 01 for 
three representative values of interaction strength. Long- 
range order can be seen in all three regimes, with C'(r) 
approaching a finite constant at large r. 

In summary, we have calculated exact properties of the 
strongly interacting 2D Fermi gas at zero temperature, by 
a combination of two AFQMC methods. The equation of 
state, contact parameter, condensation fraction and pair 
wave functions are obtained. Improved agreement is seen 
with the pressure recently measured in quasi-2D experi¬ 
ment compared to best current (approximate) theoretical 
results. Our results will provide valuable benchmarks for 
future studies and allow precise comparisons with exper¬ 
iments as the latter rapidly develop in 2D. The analytic 
forms parametrized from the accurate numerical results 
will also facilitate future local-density type of calculations 
( 4 ^ in a variety of systems relevant to experiment, includ¬ 
ing thermodynamics and out of equilibrium properties in 
the presence of a trap. The technical advances in com¬ 
putational techniques, which allowed efficient sampling 
of larger lattices with long imaginary-times and much 
smaller Monte Carlo variance than previously possible, 
can be expected to have many applications in cold atom 
systems and elsewhere. 
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FIG. 4: (Color online) Condensate fraction and pairing cor¬ 
relation functions. In the main graph, the uncertainty in the 
QMC data (from extrapolation to the TL) is estimated by 
multiple runs with different sizes and is indicated by the thick¬ 
ness of the line. Also shown are BCS results and, in the BEC 
limit, Bogoliubov results for Bose gas for reference. In the 
inset, the pairing correlation function C(r) is plotted vs. r 
for three interaction strengths (from top to bottom, the same 
parameters as in (a), (b), and (c) of Fig. |3]). The dashed lines 
are from BCS and solid lines are QMC results (error bars 
smaller than symbol size). 


Appendix A: Generalized Metropolis with force bias 


In this appendix, we describe our second approach us¬ 
ing the generalized Metropolis procedure to accelerate 
the sampling of paths in auxiliary field (AF) space. We 
introduce a dynamic force bias, analogous to what is em¬ 
ployed in the branching random walk methods in con¬ 
strained path or phase-free AFQMC [l^, in proposing 
the updates of the field values, which improves the ac¬ 
ceptance ratio and hence the MC efficiency. 

To facilitate the description of the sampling algorithm 
we first give a brief sketch of the standard path-integral 
AFQMC approach, on which more detailed descriptions 
can be found in, for example, Refs, [i^ and [l^- Ground 
state AFQMC measures the static properties by 

^ ^ (^t| exp(-/?iJ/2) 6 exp{-l3H/2) |^t) 

(V'tI exp(-/3i7) IV't) 


where the Hamiltonian H = K + V is given by Eq. ©■ 
We apply the usual Trotter-Suzuki breakup 


g-ArA g-Arkr/2g-ATCg-ATk'/2 


(A2) 


and the Hubbard-Stratonovich (HS) decomposition 


ArUni-friii _ ^ 
“ 2 




^{'IXi-ArU/2){ni^+nii-l) 


(A3) 


Xi — ±1 


1 


= 3 


Xi — ±1 

with cosh( 7 ) = exp(—Arl7/2), arriving at the form 


^-ArH 


dxp{x)B{x ), 


(A4) 


where x = The probability density 

function p(x) is uniform for the 2^® AF configurations 
under the choice of HS in Eq. (IA3I) . and the one-body 
propagator is B{x) = f([. bi{xi) 

The expression in Eq. (EH) is then re-written as a path 
integral of M = PjAr time slices. Let us consider the 
Z-th time slice, and introduce the notation 

= (V;T|B(x(“))H(x(“-i))---.B(x('+i))e“^^'^/2 
= e-'^"*/2.B(x('-i))B(x('-2))...H(x(i))|V'T), 

which are both single Slater determinant wave functions 
if we choose {iPt) to be a Slater determinant. The inte¬ 
grand of the path-integral in the denominator of Eq. (EH) 
then becomes 


W{x) = p{x){'^lJl\Y[k{x^)\’tpr) , (A5) 

i=l 


where x denotes the collection of AF at time slice 1. In 
the standard way of sampling W, one proposes to flip 
each auxiliary-field Xi one by one, and sweeps through 
X. We will update the entire configuration x (or a sub¬ 
cluster of X for very large system sizes), simultaneously. 
We define a force bias 


i'lpilA) 


(A6) 


and propose updates of the fields with the probability 
density: 


V{x) oc p(x) Y[ (A7) 

i=l 


which can be sampled directly. Detailed balance then 
leads to a Metropolis acceptance probability given by 


M(x —>• x') 


min{l, 


>V(x')iP(x)^ 

W(x)iP(x')^ 


(AS) 


Note that the probability function for proposing transi¬ 
tions does not depend on the “current” configuration of 
AF, i.e., 'P(x ^ x') = T’(x'). liV = W, all updates will 
be accepted. Because of the force bias, V approximates 
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W up to 0(\J At), leading to typically high acceptance 
ratio. 

Although we have used the discrete charge HS decom¬ 
position, the algorithm generalizes straightforwardly to 
continuous HS transformations. We comment that the 
use of the dynamic force bias in E g. (|A6I ) effectively in¬ 
troduces a background subtraction [271, [^ in the decom¬ 
position of Eq. (IA3I) . That is, if one were to employ the 
standard updating algorithms without the force bias, one 
would find Eq. dMI) much less efficient than a contin¬ 
uous charge decomposition which subtracts a constant 
background. This discrepancy in efficiency grows more 
as the system density decreases, which is especially rele¬ 
vant since the systems studied here are at the low density 
limit. (See Ref. [i^ for an analysis of the efficiency of HS 
transformations, and Ref. [1^ for discussion on how the 
dynamic force bias automatically introduces an optimal 
constant background shift.) 

Some other features of our algorithm are: 

• Since we always work in the dilute limit, the mem¬ 
ory is saved by only storing the wave function and 
calculating the Green function on the fly. We divide 
the path of M slices into y/M blocks, and only track 
one block each time. The wave function at the be¬ 
ginning of each block is stored. The largest number 
of wave functions stored in our code is ^ 2'/M. 

• The wave function is transformed between real and 
momentum space by fast Eourier transformation, so 
that all the one-body operators during projection 
are diagonal, and Green functions in different space 
are easily obtained. 

• When we only need the energy, we separate it into 
kinetic and potential energy. They are diagonal 
either in momentum or real space, where we do not 
need to calculated the whole Green function. To 
improve statistics, we measure the energy anywhere 
along the path and combine them, including the 
mixed estimator on both side. 

• The standard determinantal QMC formalism as 
sketched above turns out to have a divergence of 
the Monte Garlo variance. We discuss the variance 
problem and its solution separately elsewhere [d^. 
The solution involves the introduction of a bridge 
link, which we have implemented in the calculations 
presented here. The force bias and basic sampling 
algorithm described above remain unchanged. 

Appendix B: Extrapolation to the continuum limit 

We have described the extrapolation procedure of our 
lattice results to the continuum limit, and the subsequent 
analysis to reach the thermodynamic limit. Here we il¬ 
lustrate the finite size extrapolation in few-body systems. 

The extrapolation to the continuum limit, for a fixed 
number of particles, must be consistent and independent 


TABLE II: Data of the equation of state in Fig. [T] The in¬ 
teraction strength, given in the first column, are In(ofcF) = 
2 /-I-ln(2)/2, with y from —0.75 to 6 (in increments of 0.25 up 
to 2 / = 2, then increments of 0.5 up to 2 / = 5). 


In(afcF) 

Eqmc/Efg 

Error bar 

Ebcs/Efg 

-0.403426 

-5.512634 

0.000619 

-4.651252 

-0.153426 

-3.262997 

0.000487 

-2.427641 

0.096574 

-1.884889 

0.000325 

-1.078969 

0.346574 

-1.027841 

0.000453 

-0.260958 

0.596574 

-0.487667 

0.000335 

0.235190 

0.846574 

-0.137058 

0.000272 

0.536119 

1.096574 

0.096228 

0.000203 

0.718642 

1.346574 

0.256943 

0.000167 

0.829348 

1.596574 

0.371799 

0.000162 

0.896491 

1.846574 

0.456471 

0.000141 

0.937204 

2.096574 

0.521804 

0.000173 

0.961859 

2.346574 

0.572904 

0.000111 

0.976740 

2.846574 

0.647340 

0.000103 

0.990927 

3.346574 

0.700067 

0.000067 

0.997096 

3.846574 

0.737144 

0.000128 

0.997307 

4.346574 

0.767283 

0.000099 

0.997765 

4.846574 

0.793547 

0.000068 

1.000206 

5.346574 

0.813073 

0.000053 

1.000436 

6.346574 

0.842689 

0.000036 

1.000654 


of the type of kinetic energy dispersion. Eor a two-body 
problem on the lattice, exact results can be obtained for 
large system sizes by mapping to a one-body problem 
in the center of mass system. The results are shown in 
Fig. EJa), which fit well a 4th-order polynomial function 
in l/L. We see from the inset that the coefficient on the 
linear term is zero within numerical precision. 

We also show the finite size effect in the four-body 
problem from QMG, in Fig. EDb), reaching large lattice 
sizes. The same general behavior is seen as in the two- 
body problem. We have also studied the finite-size be¬ 
havior of the BCS solution, finding similar trends but 
with different slopes. In the many-body system, our 
QMG data are consistent with these observations as well. 
They are thus fitted with a 4th-order polynomial func¬ 
tion with a vanishing l/L coefficient, as described in the 
main text. 


Appendix C: Equation of state data 

We list the data for the equation of state in Fig. [T] 
The QMC energy data are calculated by our branching 
random walk approach with BCS trial wave functions. 
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FIG. 5: (Color online) Extrapolation of finite-size lattice to the continuum limit in few body problems. Panel (a) shows exact 
diagonalization results for the two-body problem at ln(afc_F) = 0.5, while panel (b) shows QMC solutions for the four-body 
problem at In(afcF) = 0.0. In each case, results are obtained for both the Hubbard and the quadratic dispersions. A 4th-order 
polynomial function in 1/L fits well both dispersions, and the extrapolated results in continuum limit agree well with each 
other. The insets indicate that the coefficients on 1/L are negligible in both cases. 


TABLE III: Final parameter values (full digits) in the 
parametrization [Eqs. (I3I5II 1 of the exact EOS from QMC. 


Oq 

-11.804127317953723 

a[ 

14.675499370762239 

(A 

-4.855080880566919 

ao 

-0.819842357425408 

ai 

0.1273251139440354 

02 

0.0685123559420463 

03 

-0.014505432043856327 

a4 

-0.009191602440101383 

Os 

0.004190575139056055 

Os 

-0.0006367374265820822 

ar 

3.431232818866204x10"® 

(22 

-0.060852400057644876 

V 

“3 

0.36401186693517423 

04 

-0.61531422724189 


Appendix D: Full digits for Table [T] 
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